Stationary Self-Organized Fractal Structures in an Open, 
Dissipative Electrical System 

LETTER TO THE EDITOR 



Marco Maranij |, Jayanth R. Banavar^:, Guido Caldarelli^f, Amos Maritan||, and Andrea Rinaldof 

f Dipartimento Di Ingcgneria Idraulica, Marittima e Geotecnica e Centro Internazionale di Idrologia "Dino Tonini", 
Universita di Padova, via Loredan 20, 1-35131 Padova, Italy 

^Department of Physics and Center for Materials Physics, The Pennsylvania State University, 104 Davey 
Laboratory University Park, PA 16802 USA 

% TCM Cavendish Laboratory, Madingley Road, CB3 OHE, Cambridge, UK 

|| International School for Advanced Studies, 1-34014 Grignano di Trieste, Istituto Nazionale di Fisica della Materia 
and sezione INFN di Trieste, Italy 



Abstract. We study the stationary state of a Poisson problem for a system of N perfectly conducting metal 
balls driven by electric forces to move within a medium of very low electrical conductivity onto which charges are 
sprayed from outside. When grounded at a confining boundary, the system of metal balls is experimentally known 
to self-organize into stable fractal aggregates. We simulate the dynamical conditions leading to the formation of 
such aggregated patterns and analyze the fractal properties. From our results and those obtained for steady-state 
systems that obey minimum total energy dissipation (and potential energy of the system as a whole), we suggest 
a possible dynamical rule for the emergence of scale-free structures in nature. 
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1. Introduction 

The key ingredients for the ubiquitous appearance of natural phenomena lacking a spatial and/or temporal 
scale remain a mystery. Self-organized criticalityjj], ^], ^| represents one basic attempt towards the elucidation 
of a general mechanism responsible for scale and time invariance in a variety of phenomena in diverse fields 
including physics, biology, geophysics and economics. An alternative view has been postulated on the basis 
of fractal forms occurring in natural river basins and in domain walls of random ferromagnets^, ^, ^, |?], |). 
It has been speculated that chance and necessity, defined by an imperfect strive for global optimality of 
open, dissipative systems could be ingredients for some forms of natural fractal growth. In this paper, we 
study the scale-free stationary states resulting from the dynamics of an open electric system. The steady state 
configuration corresponding to minimum total potential energy is a periodic non-fractal structure. Our findings 
suggest a rule for the emergence of fracatality that encompasses dynamical problems with a global constraint. 

We propose a model for the spontaneous aggregation of metal balls floating in a viscous liquid having a 
high dielectric constant. Electric charges sprayed on to the poorly conducting still fluid are transported to a 
grounded confining ring by the metal balls. Such aggregation has been observed]^, |l^| in cylindrical cells made 
up by a layer of castor oil (because of its high dielectric constant) within an acrylic dish confined by a layer 
of air above. The inner perimeter of the dish was equipped with a grounded metal ring. A strong potential 
between a metallic tip and this ring was fixed and charges were sprayed quasi homogeneously upon the oil 
surface. Inside the dish there were metallic bearing balls, which spontaneously arrange themselves under the 
influence of electric forces that lead to the transportation of charges to the metal ring. It is relevant to our 
computations that the reference experiment was repeated several times for different numbers of balls and, at 
stationary conditions, the aggregated structures showed fractal properties^, [jij. 



2. Description of the model 

Our model consists of N particles (the metal balls of experiment || |l(|) in a domain T with Nr (^$> N) sites 
in a 2 — d square lattice. The particles are charged (due to constant injection of charges from outside) and are 
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subject to the force due to the electric potential <p(x) which obeys the Poisson equation 

vV = -SoM) = (!) 

where Sq is the source term, ao is the specific conductivity of the medium and p is the charge density. Eq. 
(|jj) follows from charge conservation and Ohm's law || |1C[ ]. The charge density is assumed to be constant 
in space (i.e. perfectly homogenous spraying of charges) and, by a suitable rescaling of the field tp, its value 
may be assumed equal to 1 and is immaterial as long as it is not zero. The Laplacian in eq. (0) is a lattice 
discretization of the corresponding continuum operator and it is computed assuming a eight-neighbour scheme 
(four neighbours along the vertical and horizontal directions and four along the diagonals). 

Poisson's eq. (Q) is solved with the boundary condition ip{dY') = 0, where dT' includes the boundary, 
<9r, of r and all of the sites in T which are occupied by a particle connected to dT through other particles. 
The particles in dT are immobile. This assumption considerably reduces the more complex system studied in 
|l^| , and is qualitatively akin to the experimental conditions where the difference in electrical conductivity 
between the metal balls and the embedding castor oil is quite large. 

Initially the N particles are placed at random in T and dT' = dT, i.e. no particles are in contact with 
the boundary. The potential <p(x) is computed at all Nr sites by solving eq. (Q) with b.c. <p(dT) = 0. 
Particles not belonging to dV are considered immaterial to the evaluation of tp because of their assumed very 
large conductivity. The gradient V l _ l p(x) in all eight directions (/i = 1, . . . , 8), at all sites x (§t dV) is now 
computed. The largest value |V/j</?(x)| is then found and the corresponding particle is moved in the direction 
jl (in case such a move caused double occupancy of a cell the site with the next largest value of |Vy(x)| is 
selected instead; and so on). This is then repeated until a particle becomes a neighbour of the boundary dV 
and thus a part of the boundary. Whenever this occurs the potential ip(x) must be recomputed due to the 
change in boundary conditions. This process is repeated until all the particles eventually become part of the 
boundary. 

Other dynamics may be used in which all the particles are moved simultaneously in the direction of the 
maximum gradient. The stationary configurations that emerge are qualitatively similar to those obtained 
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through the procedure described here. 

It is crucial that at stationary conditions for the arbitrary i-th metal ball located at x = q,, (i.e. when 
q, = cjj = 0) the total potential energy W oc J r (pdx is at a minimum with respect of the arrangement of the 
metal balls (i.e. V qi W = 0) @, Further more, the system configuration that minimizes W also exactly 
minimizes the total energy dissipation P cx J r \ Vip\ 2 dx [To| |. The above framework relates directly to other 
contexts, where scale invariance was argued to arise as a consequence of frustrated optimality ||. 

The computation of the potential </?(x) allows the determination of the total electric resistance, say R which, 
in the experiment, had been shown to scale with the number N of metal balls as R oc . Box counting 
and cluster dimensions of the sets formed by the connected aggregates were also computed and the scaling of 
N\{1), the number of boxes of size I making up the structure, and of A^r), the number of balls inside a radius 
r was studied. 

Technically, the solution of a Poisson equation at each time step implies a considerable numerical burden 
and an efficient numerical scheme had to be employed. Finite difference approaches require the solution of 
a L 2 x I? linear system, where L is the size of the lattice. The use of the accelerated conjugate gradient 
technique |Tl|], though very efficient when a single solution is required, does not take advantage of the fact that 
two successive fields only have modest differences due to the local variation of the boundary conditions. A 
combined use of a relaxation technique and an initial direct solution proved computationally more advantageous, 
because it bypasses the slow relaxations when first solving the equations, and updates, given the modified 
boundary conditions, the configuration of ^(x) and quickly converges. Note that the improvement in efficiency 
was especially critical when exploring a substantial number of possible perturbations of a given arrangement 
obtained from the dynamics to test its optimality. 

3. Numerical experiments 

We have run several simulations of the dynamics of the system where initially the N balls are randomly 
distributed over the domain T, and have obtained consistent fractal patterns for the aggregate arrangement 
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of the grounded elements. Figure [l] shows one realization of the final potential field, in this case related to a 
square domain grounded at its center (i.e. dT consists of a single site) where N — 2048 aggregated metal balls 
are employed. We have used these particular boundary conditions, at no cost of generality, to compare the 
final structure with an exact optimal arrangement described below. The shades of gray in Figure 1 are related 
to the local value of the potential y(x). We have examined in detail the statistics of the resulting aggregates at 
different densities N/N?. Figure 2 shows sample results of finite-size scaling analyses related to the estimation 
of box-counting dimensions |T^| and the mass dimension D m [[l3| . We have also computed the scaling exponent 
£ of total resistance, R oc iV _ ^[|[ Nevertheless it is not possible to compare scaling exponents from our 
numerical simulations and observations. In fact, the dynamical model is simple and does not capture all 
the experimental details. We have made simplifying assumptions for computational convenience, e.g. infinite 
electric conductivity of the metal balls, zero mechanical inertia of the metal balls and infinite cohesive forces 
among the grounded balls. 

In the case of Figure [l] (N = 2048) for a square grid of size L = 128, the box-counting exponent is D = 1.30, 
the mass dimension D m = 1.63 (Figure |^) and the scaling exponents of resistances is £ = 2.00. In this case 
the total potential energy is W = 14.9 (in arbitrary units) and in all realizations with the same boundary 
conditions and number N the final value of W is within a few units. The statistics of the algebraic laws are 
robust. We have also tested circular and square domains grounded at the outer boundaries - one obtains 
DLA-like dendrites issuing from the center with similar results. Note that in all cases the dynamics is not 
stochastic, differently from the cases of DLA and DBMju|, |l5| . 

As one expects, the global minimum is very unlikely to be accessible by the dynamics. In fact the structure 
having minimum total potential W, which we have found exactly for regular geometries (It can be shown that 
a regular grid of N connected grounded points is the global minimum of total potential VV. The proof proceeds 
from a one dimensional situation in which the arrangement of arbitrary N grounded points within a field 
obeying ip" = 1 pinched to zero at the boundaries is a regularly spaced array of zeroes. The extension to two 
dimensions comes from separation of variables of the 2D Poisson equation and the requirement of connectivity 
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of the grounded structure), is characterized by great regularity and substantially different scaling exponents 
(Figure 3). In the case of the square domain occupied by N — 2048 connected sites of zero potential, a regular 
grounded lattice yields the optimal energy W — 5.60, and a box-counting and mass exponents D = D m = 2.00, 
outside the error bars of dynamically accessible states. 

The resulting structures, in all cases, do indeed correspond to local minima of the total potential 
energy. Hence, in the exactly known case of a regular arrangement and in analogy with other physical 
processes^, (t[ ||, imperfect optimal search procedures yield local optima showing scaling properties quite 
different from those of the true ground state. In the general case, the exact properties of the ground state 
are not known, and thus we have used refined annealing procedures [fl7| |l8| starting from the self-organized 
aggregates to search for optimal aggregation. This procedure uses a randomized approach in the exploration of 
the values of the objective function in order to avoid being trapped in local minima and it is known to ensure 
achievement of values close to the absolute minimum. The arrangements obtained through this optimal search 
also show significant departures from the features of the initial structure. 

4. Conclusions 

The numerical experiments described before seem to indicate that the type of optimization that is dynamically 
feasible in nature is rather myopic, that is, willing to accept changes only if their impact is favourable in the 
immediate rather than in the long run. Whether or not dynamical rules necessarily obey a variational principle, 
and thus whether fractal forms in nature are or are not metastable states of a dynamical search for optimality, 
our results suggest that one ingredient for the emergence of fractality is the existence of a set of stationary 
or recursive states satisfying a global constraint. In our case the constraint is the dissipation of the injected 
charge through the boundaries. It is thus possible that, in some cases, local interactions may yield a globally 
felt constraint through boundary conditions. The system studied here is quite analogous to glasses. There, the 
ground state is known to be a crystal but on rapid quenching of a liquid, the dynamics inhibits crystallization 
and a glassy structure (that is not fractal) is obtained. The key difference is that here the boundaries impose a 
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global constraint that result in the metastable states becoming scale-free over a certain range of length scales. 
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Figure 1. A realization of the potential field with N = 2048 metal balls within a 128 X 128 grid grounded at the 
central site P (i.e. ip(P) = 0). In the inset we show the planar aggregation structure of the grounded metal balls. 
The total potential energy is W = 14.0 (in arbitrary units). 

Figure 2. (a) Box-counting scaling relationships, iVi(i) oc l~ D , for a varying number of balls within a 128 X 128 
domain. Notice the consistency in the slopes; (b) the computation of the mass dimension D m requires the 
evaluation of the exponent of the scaling relationship N2(r) oc r Dm where N2(r) is the number of balls falling 
within a circle of radius r seeded at P. The computation was carried out over several realizations obtained for 
2048 metal balls on a 128 X 128 gridded domain. The departure from a power law for large r is due to finite-size 
effects. In all plots larger symbols indicate the ensemble mean while dots indicate single realizations. 

Figure 3. Potential field arising from a regular, grid-like arrangement of N = 2048 metal balls within a center- 
grounded square 128 X 128 domain. W = 5.6 (in the same units employed for the system in Figure 1) and is 
substantially lower than the values obtained by spontaneous self-organization and by imperfect optimal search 
techniques. 
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